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Combinatorial optimization problems are computationally hard in general, but they are ubiquitous 
in our modern life. A coherent Ising machine (CIM) based on a multiple-pulse degenerate optical 
parametric oscillator (DOPO) is an alternative approach to solve these problems by a specialized 
physical computing system. To evaluate its potential performance, computational experiments are 
performed on maximum cut (MAX-CUT) problems against traditional algorithms such as semidef¬ 
inite programming relaxation of Goemans-Williamson and simulated annealing by Kirkpatrick, et 
al. The numerical results empirically suggest that the almost constant computation time is required 
to obtain the reasonably accurate solutions of MAX-CUT problems on a CIM with the number of 
vertices up to 2 x 10 4 and the number of edges up to 10 8 . 


I. INTRODUCTION 

Combinatorial optimization appears in many impor¬ 
tant fields such as computer science ]1], H}, drug discovery 
and life-science [§|, and information processing technol¬ 
ogy 3 ■ One of the example of such problems is an Ising 
problem to minimize the Ising Hamiltonian, which is a 
function of a spin configuration a = (of) defined as 

'H(o) =-^JijOiOj (1 <i,j<N), (1) 

i<j 

where each spin takes binary values cq G {±1}, a real 
number symmetric matrix ,Jjj denotes a coupling con¬ 
stant, and N is the total number of spins. Despite 
its simple statement, it belongs to the non-deterministic 
polynomial-time (NP)-hard class to find the ground state 
of the Ising model on the three-dimensional lattice 
Similarly, a maximum cut (MAX-CUT) problem in the 
graph theory is to find the size of the largest cut in a 
given undirected graph. Here, a cut is a partition of the 
vertices V into two disjoint subsets (Si, S 2 } and the size 
of the cut is the total weight of edges Wij with one vertex 
i in Si and the other j in S 2 ■ The size of the cut can 
be counted by assigning the binary spin values to express 
which subset the vertex i belongs to Oi G {±1} 

ru , V" \' (1 — Oi<jj) 

c w) = 2^ w u = X w u - 2 - 
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where T~L is an Ising Hamiltonian defined in Eq. m with 
Jij = —Wij. It indicates that the MAX-CUT problem is 
equivalent to the Ising problem except for the constant 
factor. 

The MAX-CUT problem belongs to the NP-hard class 
in general, even though there are graph topologies which 


can be solved in polynomial time i5j, [7H11]. Many at¬ 
tempts have been made to approximately solve NP-hard 
MAX-CUT problems, but the probabilistically checkable 
proof (PCP) theorem states that no polynomial time al¬ 
gorithms can app roximate MAX-CUT problems better 
than 0.94118 [lJ|T3|. Currently, the approximation ra¬ 
tio of 0.87856 achieved by the Goemans-Williamson algo¬ 
rithm (GW) based on semidefinite programming (SDP) 
is the best value for performance guarantee 14]. This al¬ 
gorithm is a well-established benchmark to evaluate any 
new algorithms or computing methods. 

Besides, there exist several heuristic algorithms to 
tackle these NP-hard MAX-CUT problems. The sim¬ 
ulated annealing (SA) was designed by mimicking the 
thermal annealing procedure in metallurgy [li|. A quan¬ 
tum annealing technique was also formulated and was 
shown to have competitive performance against SA 0- 
[20| . Independently, novel algorithms which are superior 
either in its speed or its accuracy are proposed [211423]]. 

We recently proposed a novel computing system to im¬ 
plement the NP-hard Ising problems using the criticality 
of laser [23427} and degenerate optical parametric oscil¬ 
lator (DOPO) phase transition [28], |29|. The architecture 
of this machine is motivated by the principle of laser and 
DOPO in which the mode with the minimum loss rate is 
most likely to be excited first. The energy of the Ising 
Hamiltonian can be mapped onto the total loss rate of the 
laser or DOPO network. The selected oscillation mode 
in the laser or DOPO network corresponds to the ground 
state of a given Ising Hamiltonian, while the gain acces¬ 
sible to all other possible modes is depleted due to the 
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cross-gain saturation. This means that a mode with the 
lowest loss rate reaches a threshold condition first and 
clumps the gain at its loss rate, so that all the other 
inodes with higher loss rates stay at sub-threshold con¬ 
ditions. Moreover, the DOPO is in the linear superposi¬ 
tion of 0-phase state and 7r-phase state at its oscillation 
threshold [3Cj. The coupled DOPOs form quantum en¬ 
tanglement in spite of their inherent dissipative natures 
[■Til , so that some form of quantum parallel search could 
be embedded in a DOPO network. 

In this article, the validity of the CIM for MAX-CUT 
problems is tested against the representative approxima¬ 
tion algorithms. The DOPO signal pulse amplitudes in 
CIM, which are interpreted as the solution, are described 
by the c-number stochastic differential equations (CSDE) 
as presented in Section El Then we conduct numerical 
simulations for MAX-CUT problems in Section Hill with 
the number of vertices up to TV = 20000. It is, of course, 
difficult to compare the performance of the proposed sys¬ 
tem as a MAX-CUT solver with the representative ap¬ 
proximation algorithms which can be run on current dig¬ 
ital computers mainly because the unit of “clock” can¬ 
not be uniquely defined. Thus we defined the feasible 
system clock which dominates the computing process in 
CIM as mentioned later. Moreover, here we evaluated 
the computational ability under either time or accuracy 
was fixed, while a preliminary benchmark study done in 
the previous paper is focused on the performance after 
physical convergence [H]. 


II. MULTIPLE-PULSE DOPO WITH MUTUAL 
COUPLING 

A. A proposed machine 

A standard CIM based on multiple-pulse DOPO with 
all-optical mutual coupling circuits is shown in Fig. CD 
The system starts with a pulsed master laser at a wave¬ 
length of 1.56 //m. A second harmonic generation (SHG) 
crystal produces the pulse trains at a wavelength of 
0.78 fj ,m which in turn generate multiple DOPO pulses at 
a wavelength of 1.56 /im inside a fiber ring resonator. If 
the round trip time of the fiber ring resonator is properly 
adjusted to TV times the pump pulse interval, we can si¬ 
multaneously generate TV independent DOPO pulses in¬ 
side the resonator. Each of these pulses is either in 0- 
phase state or 7r-phase state at well above the oscillation 
threshold and represents an Ising spin of up or down. 

In order to implement an Ising coupling in Eq. ©, 
a part of each DOPO pulse in the fiber ring resonator is 
picked-off and fed into an optical phase sensitive amplifier 
(PSA), followed by optical delay lines with intensity and 
phase modulators. Using such TV — 1 optical delay lines, 
(arbitrary) T-th pulse can be coupled to (arbitrary) j-th 
pulse with a coupling coefficient Jy. Such an all-optical 
coupling scheme has been demonstrated for TV = 4 and 
TV = 16 CIMs Hll|. 


In Sections IIII Dl and IIIIE1 we assume a CIM with a 
fiber length of 2 km (or cavity round trip time of 10 ps) 
and pulse spacing of 10 cm (or pulse repetition frequency 
of 2 GHz), thus 2 x 10 4 independent DOPO pulses can be 
prepared for computation. The system clock frequency 
for the CIM should be defined by the cavity circulation 
frequency (inverse of cavity round trip time). One clock 
cycle (round trip) includes every elements of computa¬ 
tion, such as parametric amplification, out-coupling port, 
and coherent feedback. Thus the clock frequency of the 
CIM assumed for the present benchmark study is 100 
kHz since the round trip time of 2 km fiber ring is 10 /rs. 
We fixed this system clock frequency, just like any dig¬ 
ital computer has a fixed clock frequency and chose the 
appropriate pulse interval to pack the desired number of 
pulses in the 2km fiber. 


B. c-number stochastic differential equations for 
multiple-pulse DOPO with mutual coupling 

The in-phase and quadrature-phase amplitudes of a 
single isolated DOPO pulse obey the following c-number 
stochastic differential equations (CSDE) [33]: 

dc = (—1 + p — c 2 — s 2 )cdt + -j- \ic 2 +s 2 + ^-d,Wi (3) 

A s V 2 

ds = (—1 — p — c 2 — s 2 )s dt + -j -1 /c 2 + s 2 + ^(iW 2 .( 4 ) 

A. s V z 

The above CSDE are derived by expanding the DOPO 
field density operator with the truncated Wigner dis¬ 
tribution functions. An alternative approach is to use 
two coherent states in the generalized (off-diagonal) 
P(a s , ^-representation for the field density matrix JUj. 
The two approaches by the truncated Wigner function 
and the generalised P-representation are equivalent for 
highly dissipative systems such as ours. The pump 
field is adiabatically eliminated in © and © by as¬ 
suming that the pump photon decay rate y p is much 
larger than the signal photon decay rate 7 S . The term 
A s = is the DOPO field amplitude at a 

normalized pump rate p = F p /F t h = 2, and k is the sec¬ 
ond order nonlinear coefficient associated with the de¬ 
generate optical parametric amplification. The variable 
t = 7 s r /2 is a normalized time, while r is a real time 
in seconds. The term F p is the pump held amplitude 
and Fth = 7 s 7 p/ 4 k is the threshold pump held ampli¬ 
tude. Finally, dW\ and dW 2 are two independent Gaus¬ 
sian noise processes that represent the incident vacuum 
fluctuations from the open port of the output coupler and 
the pump held huctuation for in-phase and quadrature- 
phase components, respectively. The vacuum huctuation 
of the signal channel contributes to the 1/2 term and the 
quantum noise of the pump held contributes to c 2 + s 2 
in the square-root bracket in © and ©. 

When the T-th signal pulse is incident upon the out¬ 
put coupler, the output-coupled held and remaining held 
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FIG. 1. A coherent Ising machine based on the time-division multiplexed DOPO with mutual coupling implemented by optical 
delay lines. A part of each pulse is picked off from the main cavity by the output coupler followed by an optical phase sensitive 
amplifier (PSA) which amplifies the in-phase amplitude Ci of each DOPO pulse. The feedback pulses, which are produced 
by combining the outputs from N — 1 intensity and phase modulators, are injected back to the main cavity by the injection 
coupler. 


inside a cavity are written as 


Ci ,out = VTcj -VI -T-J- 
A s 

(5) 

Ci, cavity = Vl-Ta+VT^-, 

( 6 ) 


where T is the power transmission coefficient of the out¬ 
put coupler and fa is the incident vacuum fluctuation 
from the open port of the coupler. The out-coupled field 
and the signal amplitude after PSA can be 


~ Ci out 

Ci = r— — Q 

VT 


1 -T fa 
T V 


(7) 


From these out-coupled pulse stream, the intensity and 
phase modulators placed in the N — 1 delay lines pro¬ 
duce the mutual coupling pulse ]Th fa,jCj- which is actu¬ 
ally added to the i-th signal pulse by an injection coupler. 
Here, faj is the effective coupling coefficient from the j-th 
pulse to the i-th pulse, determined by the transmission 
coefficient VT' of the injection coupler. In the highly dis¬ 
sipative limit of a mutual coupling circuit, such as in our 
scheme, we can use the CSDE supplemented with the 
noisy coupling term. Since the transmission coefficient 
VT' of the injection coupler should be much smaller than 
one, we do not need to consider any additional noise in 
the injected feedback pulse. The CSDE ([3]) can be now 
rewritten to include the mutual coupling terms 


da = [(-1 +P-C-- Si)a + ^2 dt 

3 

+ J~ s \l C i +S i + \ dWi - ( 8 ) 

The summation in Eq. © represents the quantum 
measurement-feedback term including the measurement 
error given by Eq. 0 . The vacuum fluctuation coupled 
to the i-th pulse in the output coupler is already taken 
into account in the last term of right-hand side of Eq. 


© together with the pump noise. We conducted the nu¬ 
merical simulation of the coupled CSDE © to evaluate 
the performance of the CIM. 


III. BENCHMARK STUDIES ON MAX-CUT 
PROBLEMS 

A. MAX-CUT problems on cubic graphs 

The MAX-CUT problem on cubic graphs, in which 
each vertex has exactly three edges, is called MAX- 
CUT-3 problem and also belongs to NP-hard class [35 j. 
The smallest simple MAX-CUT-3 problem is defined on 
the complete graph K 4 with four vertices and six edges 
with identical weight Jij = — 1 , where anti-ferromagnetic 
couplings have frustration so that the ground states 
are highly degenerate. The solution to this problem 
are the set of two-by-two cuts, which contains six de¬ 
generate ground states of the Ising Hamiltonian, i.e., 

{ltm>, mu), imt), \xm ), imt), iutt». 

Figure [5] shows the time evolution of a (* = 1, ■ ■ •, 4) 
when p = 1.1 and £ = — 0.1. A correct solution spon¬ 
taneously emerges after several tens of round trips. The 
statistics of obtaining different states against 1000 ses¬ 
sions of such a numerical simulation are shown in Fig. [3l 
in which six degenerate ground states appear with almost 
equal probabilities with no errors found. 


B. Many-body interaction problem 

If the interaction is not a standard two-body Ising in¬ 
teraction type but rather a four-body interaction such 
as 

H = — Jl2340’l0'2CT3tT4, (9) 

where J1234 € R, the coupled field into the i-th pulse is 
no longer given by but b y £cjC k ci (j,k,l fa i). 
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FIG. 2. Normalized DOPO signal amplitudes as a function of 
normalized time (in unit of cavity round trip numbers) for a 
N = 4 simple MAX-CUT-3 problem. Each color corresponds 
to the four different DOPOs indexed with i = 1,..., 4. Small 
window is enlarged to indicate the status of signal amplitude 
inside a cavity at three components (as in Fig. [TJ; A: OPA 
gain medium (PPLN waveguide), B: out-coupler, and C: in¬ 
jection coupler of the mutual coupling pulse. The two flat 
regions between B and C and between C and A are the pas¬ 
sive propagation in a fiber. 
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FIG. 3. Distribution of output spin configurations in 1000 
trials of numerical simulations against a simple MAX-CUT-3 
problem of graph order N = 4. All trials were successful to 
find one of the six degenerate ground states. 
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FIG. 4. Normalized DOPO pulse amplitudes d (i = 1,... ,4) 
under the interaction between four-body Ising coupling ex¬ 
pressed by Eq. ©• 
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FIG. 5. Distribution of output spin configurations in 1000 
trials of numerical simulation against a four-body Ising model 
of N = 4. All trials were successful to find one of the eight 
degenerate ground states. 


In this case, the CSDE m can be rewritten to include 
the four-body coupling term 

dci = [(-1 + p - cf - sl)ci + ticket] dt 

+ J~ s \J c i +s i + \ dWi - ( 10 ) 

When the four-body coupling coefficient J1234 
is —1 (multi-body anti-ferromagnetic coupling), 
there are eight degenerate ground states, i.e., 
Itttl), Ittlt), Itltt), 14-ttt) and their inverse spin 
configurations. Figure [I] shows the time evolution of Ci 
[i = 1,... ,4) when p = 1.1 and £ = —0.1. One of the 
eight degenerate ground states emerges spontaneously 
after several tens of round trips. The statistics of 
observing different states in 1000 independent sessions of 
the numerical simulation of Eq. dll are shown in Fig. 


C. Algorithm description 

In this subsection, we will review the four representa¬ 
tive approximation algorithms for MAX-CUT problems. 

The Goemans-Williamson algorithm (GW) based on 
SDP has a 0.87856-performance guarantee for NP-hard 
MAX-CUT problems (I3 |. It achieves the optimal ap¬ 
proximation ratio for MAX-CUT problems under the 
assumptions of P ^ NP and the unique games conjec¬ 
ture [3j|. The SDP relaxation of the original MAX-CUT 
problem is a vector-valued optimization problem as max¬ 
imizing 3 Y^i<j w iA^-Vi-Vj), Vi S where is 

a unit sphere in R fc and k < N (or number of ver¬ 

tices). There exist polynomial time algorithms to find the 
optimal solution of this relaxation problem (with error 
e > 0), and its value is commonly called the SDP upper 
bound. A final solution to the original MAX-CUT prob- 
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lem is obtained by projecting the solution vector sets to 
randomly chosen one-dimensional Euclidean spaces (i.e., 
dividing the sphere by random hyperplanes). 

There are three types of computational complexities 
of the best-known algorithms for solving the SDP relax¬ 
ation problem. If a graph with N vertices and m edges is 
regular, the SDP problem can be approximately solved in 
almost linear time as 0(m) = 0(m( log N) 2 e~ 4 ) using the 
matrix multiplicative weights method [37}, where e rep¬ 
resents the accuracy of the obtained solution. However, 
slower algorithms are required for general graphs. If the 
edge weights of the graph are all non-negative, the fastest 
algorithm runs in 0(Nm ) = 0(Nm(\ogN) 2 £~ 3 ) time 
based on the Lagrangian relaxation method [38]. For 
graphs with both positive and negative edge weights, the 
SDP problem is commonly solved using the interior-point 
method, which scales as 0(N 35 ) = 0(N 35 log(l/e)) 
[39} . Besides, low rank formulation of SDP is effec¬ 
tive when the graph is sparse [4(il - [42[ . In our compu¬ 
tational experiments, the COPL_SDP based on the inte¬ 
rior point method was used as a solver for MAX-CUT 
problems [43]. The SDP upper bound C/sdp and the so¬ 
lution Cqw were obtained using the following parame¬ 
ters: interior point method was used until the relative 
gap r gap = 1 - Pobj/D 0 bj reached 10” 3 , where P Q bj and 
P obj are the objective functions of the primal and dual of 
the SDP problem, respectively [i3| . Random projection 
onto one-dimensional space was executed N times. 

For many practical applications, heuristic algorithms 
are more convenient to use, since the GW algorithm gen¬ 
erally requires long computation time 0(N 3 - 5 ). Metropo¬ 
lis et al. introduced a simple algorithm that can be used 
to provide an efficient simulation of a collection of atoms 
in equilibrium at a given temperature [45 ] . Kirkpatrick 
et al. applied the algorithm to optimization problems 
by replacing the energy of the atomic system to the cost 
function of optimization problems and using spin config¬ 
urations <7, which is called the simulated annealing algo¬ 
rithm (SA) [lj| . In each step of this algorithm, a system 
is given with a random spin flip and the resulting change 
AE in the energy is computed. If AE < 0, the spin 
flip is always accepted, and the configuration with the 
flipped spin is used as the starting point of the next step. 
If AE > 0, the spin is treated probabilistically, i.e., the 
probability that the new configuration will be accepted 
is P(AE ) = exp(—AE/kBT) with a control parameter 
of system temperature T. This choice of P(AE) results 
in the system evolving into an equilibrium Boltzmann 
distribution. Repeating this procedure, with the tem¬ 
perature T gradually lowered to zero in sufficiently long 
time, leads spins a to convergence to the lowest energy 
state. In practical case, with the finite time, the anneal¬ 
ing schedule affects the quality of output values. Here in 
our numerical simulations, the temperature was lowered 
according to the logarithmic function [46}. Note that 1 
Monte Carlo step corresponds to N trials of spin flip. 

Sahni and Gonzalez constructed a greedy algorithm for 
MAX-CUT problems, which has 1/2-performance guar¬ 


antee [2l|, and SG3 is a modified version of it [22}. In 
this algorithm, nodes V are divided into two disjoint 
subsets {Si, £2} sequentially. For each iterative pro¬ 
cess, the node with the maximum score is selected, and 
it is put into either set Si or £2 so as to earn larger 
cuts. Here, the score function of SG3 is defined as 
Xi = | Ejes, W H ~ E je s 2 w iJ | (i = U • • ■, N). It stops 
when all the edges are evaluated to calculate the score 
function, thus SG3 scales as 0{m). 

The power of breakout local search (BLS) appears in 
the benchmark result for G-set graphs [23}. It updated 
almost half of the best solutions in G-set with the special¬ 
ized data structure for sorting and dedicated procedure 
to escape from local minima. The algorithm is combi¬ 
nation of steepest descent and forced spin flipping: after 
being trapped by a local minima as a result of steep¬ 
est descent procedure, three types of forced spin flipping 
(single, pair, and random) are probabilistically executed 
according to the vertex influence list (i.e., which vertex 
will increase the number of cut most when it’s flipped) 
on each subset of partition. 

These algorithms are coded in C/C-l—|- and run on a 
single thread of a single core on a Linux machine with 
two 6-core Intel Xeon X5650 (2.67 GHz) processors and 
94 GB RAM. The CIM is simulated based on the coupled 
CSDE {8} on the same machine. Note that the computa¬ 
tion time of CIM does not mean the simulation time on 
the Linux machine but corresponds to the actual evolu¬ 
tion time of a physical CIM. 


D. Computational accuracy on G-set instances 

The performance of a CIM with DOPO network 
was tested on the NP-hard MAX-CUT problems on 
sparse graphs, so-called G-set [47|. These test instances 
were randomly constructed using a machine-independent 
graph generator written by G. Rinaldi, with the number 
of vertices ranging from 800 to 20000, edge density from 
0 .02% to 6%, and topology from random, almost planar, 
to toroidal. 

The output cut values of running the CIM, SA, GW, 
and the best known solutions so far we could find [23;, [401 
for some of G-set graphs are summarized in Tabled] The 
results for CIM are obtained in 50 ms, which correspond 
to the performance of an experimental system after 5000 
DOPO cavity round trips. The best result and ensemble 
average value for 100 trials are shown. Here, the param¬ 
eters are set to be p = 1.6, £ = —0.06 and the coupling 
constant /y/(k ) is normalized by the square 

root of the graph average degree (k). The hysteretic op¬ 
timization method, in which the swinging and decaying 
Zeeman term that flips the signal amplitude (spin) back 
and forth, is implemented four times after 10 ms initial 
free evolution [49|. Each hysteretic optimization takes 10 
ms so that the total search takes 50 ms. The result of 
SA is also obtained in 50 ms for each graph. For GW, 
the computation time ranged between 2.3 s and 1.1 x 10 5 
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s, depending on N. The best outputs of the CIM were 
1.62 ± 0.58 % better than GW but 0.38 ± 0.40 % worse 
than SA, and CIM found better cut against GW except 
for a toroidal graph (g50) and a disconnected random 
graph (g70). 

As the size of optimization problems increases, the av¬ 
erage accuracy is important for practical applications. 
Table U shows that for all G-set graphs, the average ac¬ 
curacy in 100 trials is 0.94148 to the SDP upper bound 
[/sdp, i.e., the CIM can find a cut value larger than 
0.94148 of the optimal value for the MAX-CUT problems 
on average, whereas the average accuracy of the GW is 
0.93025 and that of the SA is 0.94692. Note that Usdp 
is always greater than or equal to the optimal value for 
each MAX-CUT problem. 


E. Computation time on MAX-CUT problems 

In the previous section, the running time of CIM and 
SA are fixed to estimate the computational accuracy. 
These two algorithms explored the solutions as good as 
possible in 50 ms. Although, if we finish the computa¬ 
tion at a certain accuracy, more reasonable computation 
time can be defined. Here, the GW solution was used 
as the mark of sufficient accuracy because it ensures the 
87.856% of the ground states. The CIM and SA then 
competed the computation time to reach the same values 
obtained by GW. The time and temperature scheduling 
parameters of the SA were set as follows: Inverse temper¬ 
ature increased with logarithmic function. The number 
of spin flipping was optimized to be 10 ; times for some 
l £ N, which requires the minimum computation time to 
achieve the same accuracy as with the GW. 

Computational experiments were conducted on fully 
connected complete graphs, denoted by AW, where the 
number of vertices N ranging from 40 to 20000 and the 
edges are randomly weighted ±1. Figure [6] shows the 
Ising energy in Eq. m as a function of running time. 
Both CIM and SA run stochastically due to quantum 
and thermal noise, respectively, the ensemble average of 
energies are calculated as follows: For the CIM, the en¬ 
ergy of all 100 runs was averaged at each round trip. 
For the SA, the averaged energy was calculated at each 
point on the time axis with an interpolated value from 
real time sampling. The parameters for CIM are chosen 
to be p = 0.2, £ = -0.03 (for N = 800), and £ = -0.003 
(for N = 4000). In Fig. El (a), where N = 800, the GW 
achieved an energy equal to —15624 in 22.96 s, while the 
CIM and SA reached the same energy in 6.1 x 10~ 4 s and 
in 0.671 s. In Fig. [6] (b), the GW achieved an energy of 
— 168160 in 6646.25 s, while the CIM reached the same 
energy in 5.5 x 10 -4 s and the SA did so in 10.1 s. 

Note that this result of SA and GW comes from a spe¬ 
cific computer configuration as mentioned in Sec. IIII Cl 
There is room for an improvement in the computation 
time in constant factor due to cases like using faster 
CPUs or parallelized codes. Similarly, the computation 


time of CIM also depends on the system configuration 
and can be made faster when we use the higher clock 
frequency. In this sense, the ratios between time of CIM 
and that of the other algorithms are arbitrary. Thus we 
should study the computation time scaling as a function 
of the problem size. 

Figure [7] (a) shows the computation time versus prob¬ 
lem size (number of vertices). The computation time is 
defined as the CPU time to solve a given MAX-CUT 
problem in complete graph for GW; as the CPU time to 
reach the same accuracy as GW for SA, SG3, and BLS; 
and as the time estimated by the (number of round trips) 
x (cavity round trip time) to obtain the same accuracy 
as GW for CIM. The preparation time needed to input 
Jij into the computing system, i.e., the graph I/O time, 
is not included. For complete graphs of N < 20000, the 
CIM exhibits a problem-size independent computation 
time of less than 10 ~ 3 s if we assume the fixed cavity 
circulation frequency of 100 kHz and pulse interval of 
10 cm. This means the target accuracy is obtained in 
the constant number of round trips. It indicate that the 
computation time of CIM is determined by the turn-on 
delay time of the DOPO network oscillation, which in 
turn depends on the round trip time and the pump rate. 

In Figurc[7|(b), the computation time of the CIM with 
different system clock frequency and pulse spacing are 
shown (see the Sec. Ill Al for the definition). Since the 
solutions are obtained in a constant number of cavity 
round trips, the computation time is pulse spacing in¬ 
dependent but linearly depends on the clock frequency, 
i.e., cavity circulation frequency. The number of pulses 
accommodated in the fiber can be changed to vary the 
pulse spacing under the fixed clock frequency. On the 
other hand, when the pulse spacing is fixed and the fiber 
length is varied, the maximum number of pulses should 
be increased in proportional to the fiber length. 

Going back to the Figure [7] (a), the time complexity 
0(N 3 5 ) for the GW is dominated by the interior-point 
method in the Goemans-Williamson algorithm. The SA 
seems to scale in 0(N 2 ), which indicates that it requires 
the number of spin flips to be proportional to N (i.e., 
constant Monte Carlo steps) to achieve the optimal per¬ 
formance. Each spin flip costs a computation time pro¬ 
portional to the degree fc,:, where ki is equal to N — 1 
for all i £ V in the case of complete graphs. Thus, 
the computation time scales as 0(N{k}) = 0(N 2 ) for 
the SA in the complete graphs. Note that CIM and SA 
didn’t always reach the energy obtained by GW for the 
graph of N = 40, half of the 100 runs of stochastic al¬ 
gorithms were post-selected to reach that value. SG3 
scales as 0(m) = 0(N 2 ) in Fig. [7] (a), but the values for 
N = 40,160,800 are not shown because it didn’t reach 
the accuracy reached by the GW solution. BLS exhibits 
competitive performance against SA. Besides, the DOPO 
amplitudes in CIM evolve as in Fig. [8] (a) when N = 800. 
The distribution of computation time for 100 randomly 
weighted complete graphs of N = 800 is also shown in 
Fig. [5](b). 
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TABLE I. Performance of the coherent Ising machine, simulated annealing and Goemans-Williamson SDP algorithm in solving 
the MAX-CUT problems on sparse G-set graphs. is the number N of vertices in the graph, is the number m of edges, 
Us dp is the optimal solution to the semidefinite relaxation of the MAX-CUT problem when the duality gap reaches to 1CU 3 , 
and Cbest is the best known result so far we could find. Cgw is the best solution obtained by N projections after SDP. Csa 
and (Csa) are the best and average values obtained by SA in 100 trials of 50 ms. Ccim and (Ccim) are the best and average 
values in CIM in 100 runs of 50 ms (= 5000 DOPO cavity round trips), respectively. To make comparisons with each other, 
every cut value C generated from each algorithm is normalized according to (C + E neg )/ (Usdp + E nes ), where E ne g > 0 is the 
number of negative edges. In the bottom of this table, the average and worst values of all 71 G-set graphs are shown. 


Graph 

#v 

#E 

UsDP 

Cbest 

Cgw 

Csa 

(Csa) 

Ccim 

(Ccim) 

gl 

800 

19176 

12083 

0.9620 

0.9457 

0.9620 

0.9597 

0.9614 

0.9570 

g6 

800 

19176 

2656 

0.9607 

0.9448 

0.9606 

0.9592 

0.9601 

0.9559 

gll 

800 

1600 

629 

0.9540 

0.9327 

0.9526 

0.9478 

0.9455 

0.9370 

g!4 

800 

4694 

3191 

0.9602 

0.9336 

0.9580 

0.9544 

0.9514 

0.9472 

g!8 

800 

4694 

1166 

0.9500 

0.9282 

0.9492 

0.9439 

0.9434 

0.9372 

g22 

2000 

19990 

14136 

0.9450 

0.9191 

0.9445 

0.9409 

0.9405 

0.9361 

g27 

2000 

19990 

4141 

0.9435 

0.9174 

0.9422 

0.9400 

0.9390 

0.9356 

g32 

2000 

4000 

1567 

0.9559 

0.9272 

0.9508 

0.9478 

0.9424 

0.9384 

g35 

2000 

11778 

8014 

0.9588 

0.9292 

0.9551 

0.9523 

0.9471 

0.9438 

g39 

2000 

11778 

2877 

0.9464 

0.9226 

0.9431 

0.9399 

0.9364 

0.9318 

g43 

1000 

9990 

7032 

0.9471 

0.9292 

0.9471 

0.9439 

0.9458 

0.9396 

g48 

3000 

6000 

6000 

1.0000 

1.0000 

1.0000 

0.9919 

1.0000 

0.9747 

g51 

1000 

5909 

4006 

0.9606 

0.9333 

0.9583 

0.9544 

0.9506 

0.9468 

g55 

5000 

12498 

11039 

0.9325 

0.9006 

0.9264 

0.9215 

0.9193 

0.9160 

g57 

5000 

10000 

3885 

0.9561 

0.9237 

0.9496 

0.9473 

0.9419 

0.9384 

g59 

5000 

29570 

7312 

0.9440 

0.9148 

0.9376 

0.9356 

0.9308 

0.9288 

g60 

7000 

17148 

15222 

0.9313 

0.8989 

0.9231 

0.9201 

0.9191 

0.9152 

g64 

7000 

41459 

10466 

0.9440 

0.9143 

0.9347 

0.9324 

0.9320 

0.9299 

g67 

10000 

20000 

7744 

0.9554 

0.9215 

0.9480 

0.9459 

0.9411 

0.9388 

g70 

10000 

9999 

9863 

0.9674 

0.9633 

0.9523 

0.9479 

0.9515 

0.9482 

g81 

20000 

40000 

15656 

0.9551 

0.9195 

0.9187 

0.9125 

0.9393 

0.9376 

Average 




0.9540 

0.9303 

0.9502 

0.9469 

0.9464 

0.9415 

Worst 




0.9313 

0.8989 

0.9187 

0.9125 

0.9185 

0.9149 


(a) (b) 




FIG. 6. Performance comparison of CIM, SA, and GW in solving complete graphs (a) ATsoo and (b) A' 4000 , where each edge 
was randomly weighted ±1. Each time, bundle of curves depicted average energy (solid line) ± standard deviations (dashed 
line) in 100 runs. Dotted line is the accuracy obtained by GW algorithm, whose computation time is shown with dot. The 
number of spin flip in SA algorithm were 10 5 for Kgoo and 10 6 for A' 4000 , respectively, to optimize the computation time. SA 
and GW are running on 2.67 GHz Intel Xeon without parallelization and CIM has the cavity circulation frequency of 100 kHz. 















(b) 



Problem size N 



FIG. 7. (a) Computation time of coherent Ising machine (CIM) empirically scales as 0(1), while that of simulated annealing 
(SA), and Goemans-Williamson SDP (GW) fitted well to lines indicating 0(N 2 ) and 0(A 3 ' 5 ), respectively. Computation time 
of SG3 and BLS also fit to 0(N 2 ) with a difference of constant factor. The computation time of CIM, SA, SG3, and BLS is 
defined as the time to reach the same accuracy achieved by GW (without I/O time). Data points of CIM, SA, SG3, and BLS 
were calculated by averaging 100 runs. Only for N = 800 case, 100 randomly ±1 weighted complete graphs are generated. The 
results for N = 800 show the average computation time for 100 graphs with error bars indicating standard deviations (which 
are distributed as in Fig. |8](b)). Note that the computation time of CIM is evaluated by (number of round trips) x (10 /rs) as 
in Sec. Ill Al (b) The computation time of CIM when the DOPO cavity circulation frequency {10 kHz, 100 kHz, 1 MHz} and 
pulse interval {10 cm (solid line), 1 cm (dashed line), 1 mm (dotted line)} are changed. Computation time of CIM scales as 
O(N) if we vary the fiber length proportional to N with fixed pulse interval. 

(a) (b) 


0.5 
0.4 
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0.0 0.5 1.0 1.5 2.0 2.5 3.0 
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FIG. 8. Detailed data for N = 800 randomly weighted complete graphs (appeared in Fig. [6] (a)), (a) Time evolution of DOPO 
amplitudes of the graph of N = 800. Only 100 pulses out of N = 800 signal pulses are shown. (b)Histogram of computation 
times for solving 100 randomly weighted complete graphs of size N = 800 by the CIM, SA, and GW. 



Computation time (s) 



Computation time for the random graphs in G-set in¬ 
stances is also studied. Here the subset of graphs in which 
MAX-CUT problems can be solved in polynomial time 
(i.e., planar graphs, weakly bipartite graphs, positive 
weighted graphs without a long odd cycle, and graphs 
with integer edge weight bounded by N and fixed genus) 
are excluded. The execution time of CIM is evaluated un¬ 
der the machine spec described in S ec. Ill Al with p = 0.2, 
£ = —0.06, and £y = £wij /y/{k) . Again, the compu¬ 
tation time of SA and CIM is the actual time (without 
graph file I/O) to obtain the same accuracy of solution as 
GW. Figure 0 shows the computation time as functions 
of the problem size N. The computational cost of interior 
point method dominates the GW algorithm. (Note that 
G-set contains graphs with both positive and negative 


edge weights so that we must use the slowest interior 
point method.) Then the computation time is almost 
constant for both SA and CIM. The computation time of 
SA with constant Monte Carlo step is expected to scale 
0(N(k)) ~ 0(1) (here for the random graphs in G-set, 
(k) ~ 0(7V -1 - 09 )). The computation time of CIM here 
is governed by a turn-on delay time of the DOPO net¬ 
work to reach a steady state oscillation condition, which 
is constant for varying values of N as mentioned above 
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FIG. 9. Computation time of coherent Ising machine, sim¬ 
ulated annealing algorithm, and Goemans-Williamson SDP 
algorithm on random graphs in G-set instances. The compu¬ 
tation time for CIM and SA is defined as the time to reach 
the same accuracy achieved by GW (without I/O time). Note 
that the computation time of CIM is evaluated by (number 
of round trips) x (10 /as) as in Sec. Ill Al 


IV. SUMMARY AND DISCUSSION 


The potential for solving NP-hard problems using a 
CIM was numerically studied by conducting computa¬ 
tional experiments using the MAX-CUT problems on 
sparse G-set graphs and fully connected complete graphs 
of order up to 2 x 10 4 . With the normalized pump rate 
and coupling coefficient p = 0.2 and £ = —0.06, the CIM 
achieved a good approximation rate of 0.94148 on aver¬ 


age and found better cut compared to the GW for 69 
out of 71 graphs in G-set. The computation time for this 
sparse graph set, including few sessions of hysteretic op¬ 
timization, is estimated as 50 ms. The time scaling was 
also tested on complete graphs of number of vertices up 
to 2 x 10 4 and number of edges up to 10 s . The results 
imply that CIM achieves empirically constant time scal¬ 
ing in a fixed system clock frequency, i.e., the fixed cavity 
circulation frequency (fiber length), while SA, SG3, and 
BLS scale as 0(N 2 ) and GW scales as 0(N 3 - 5 ). Those 
results suggest that CIM may find applications in high¬ 
speed computation for various combinatorial optimiza¬ 
tion problems, in particular for temporal networks. 

The present simulation results do not mean that the 
CIM can get a reasonably accurate solution by a constant 
time for arbitrary large problem size. As mentioned al¬ 
ready, in the CIM based on a fiber ring resonator, the 
number of DOPO pulses is determined by the the fiber 
length and the pulse spacing. In order to implement 
2 x 10 5 and 2 x 10 6 7 8 9 10 11 12 DOPO pulses in the 20 km fiber 
ring cavity, we must use a pulse repetition frequency to 
2 GHz and 20 GHz, respectively. This is a challenge for 
both optical components and electronic components of 
CIM, but certainly within a reach in current technolo¬ 
gies. 
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